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Abstract 

We consider solid state phase transformations that are caused by free energy densities with domains of 
non-convexity in strain-composition space; we refer to the non-convex domains as mechano-chemical spinodals. 
The non-convexity with respect to composition and strain causes segregation into phases with different crystal 
structures. We work on an existing model that couples the classical Cahn-Hilliard model with Toupin’s theory of 
gradient elasticity at finite strains. Both systems are represented by fourth-order, nonlinear, partial differential 
equations. The goal of this work is to develop unconditionally stable, second-order accurate time-integration 
schemes, motivated by the need to carry out large scale computations of dynamically evolving microstructures 
in three dimensions. We also introduce reduced formulations naturally derived from these proposed schemes 
for faster computations that are still second-order accurate. Although our method is developed and analyzed 
here for a specific class of mechano-chemical problems, one can readily apply the same method to develop 
unconditionally stable, second-order accurate schemes for any problems for which free energy density functions 
are multivariate polynomials of solution components and component gradients. Apart from an analysis and 
construction of methods, we present a suite of numerical results that demonstrate the schemes in action. 


1 Introduction 

Many multicomponent solids undergo phase-transformations in which a diffusional redistribution of their dif¬ 
ferent components is coupled with a structural change of the crystallographic unit cell. One example is the 
phase-transformation of yttria-stabilized zirconia Zti- x Y x 02- x /2 from cubic at high-Y composition to tetrago¬ 
nal at low-Y if quenched to a low temperature. Another is observed in the spinel structures of Li-ion electrodes 
when cubic LiMn 04 transforms into tetragonal Li 2 Mn 204 upon discharging to low voltages. Pure Zr02 and 
TaS 2 are other materials that may be susceptible to such mechano-chemical phase transformations. 

The underlying phenomenology can be described by a free energy density function, which is non-convex in 
its mechanical (strain) and chemical (composition) arguments. Cahn and Hilliard [3] famously introduced the 
chemical spinodal as the domain in composition space where the free energy density is non-convex, and varies 
smoothly between local minima that correspond to distinct phases. This notion has recently been extended 
to the mechano-chemical spinodal, defined as the domain in strain-composition space where the Hessian of the 
(sufficiently smooth) free energy density function has non-positive eigenvalues [21]; see Fig. [l] If the state of a 
material point lies within the mechano-chemical spinodal, and specifically if the free energy density is non-convex 
with respect to composition, the solid will undergo diffusional segregation. Here we concern ourselves with cases 
in which the two resulting phases have cubic and tetragonal crystal structures, respectively. The cubic structure 
corresponds to a minimum of the free energy density function in strain-composition space. Furthermore, we 
adopt the undistorted cubic structure as the reference state for strain. Then, the tetragonal lattice is naturally 
obtained by the strain relative to the cubic structure. The symmetries of the cubic lattice split into three 
identical sub-groups, each of which corresponds to a tetragonal lattice oriented along one of the cubic crystal 
axes. These three tetragonal variants, however, correspond to different strains relative to the reference cubic 
lattice. The free energy density function admits three additional minima, each corresponding to the strain 
that transforms the reference cubic lattice into one of the tetragonal variants. However, the compositions at 
these additional minima are identical. As the states of material points traverse such a multi-welled free energy 
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density surface in strain-composition space the solid develops a microstructure. The mechano-chemical spinodal 
is regarded as a domain of instability since small fluctuations in strain and composition tend to grow as the 
state evolves towards one of the wells in strain-composition space. Stress softening and “uphill” diffusion result, 
respectively. 

A non-convex free energy density function can lead to microstructure as explained above. However, a 
mathematical model restricted to the above phenomenology leads to ill-posed partial differential equations 
(PDEs) for elasticity and transport characterized by spurious mesh dependence. The volume fractions of the 
cubic and tetragonal phases would be set by initial and boundary conditions on the transport problem, and 
the microstructural pattern of tetragonal variants would depend on mechanical boundary conditions. However, 
in the absence of intrinsic length scales in the mathematical model, the extent and thickness of interfaces 
between phases and variants would be determined by the mesh. There is a well-understood physical aspect 
to this argument, also: The model promotes free energy minimizing microstructures, but incurs no penalty 
for the strain and composition gradients as these fields fluctuate between one phase or variant and the next. 
Arbitrarily fine energy minimizing microstructures are therefore admissible—an essentially unphysical result. 
Mathematical well-posedness and physical realism are restored by extending the free energy density function to 
include a dependence on gradient fields of strain and composition. The corresponding free energy coefficients 
introduce intrinsic length scales, and the gradient energies distinguish between microstructures of differing 
fineness, penalizing those that vary rapidly. 

Gradient free energies in classical settings lead to the Cahn-Hilliard equation for mass transport [3], and 
variants of strain gradient elasticity that represent size effects. Because the transformation strains between the 
cubic/tetragonal phases and between the tetragonal variants are of finite magnitude in many material systems, 
we are led to Toupin’s theory of nonlinear (finite strain) gradient elasticity [23]. The Cahn-Hilliard equation 
and Toupin’s strain gradient elasticity present fourth-order spatial derivatives in primal strong form, with 
corresponding weak forms carrying second-order derivatives. Carrying out large scale computations of these 
dynamically evolving microstructures in three dimensions is challenging; efficient and accurate time-integration 
algorithms are demanded. Rudraraju and co-workers used the Backward Euler algorithm m to solve the same 
problem that we consider here. However, in that work the authors concerned themselves with introducing the 
notion of the mechano-chemical spinodal, and exploring the associated physics; not with a development/an 
analysis of accurate schemes, which is the goal of the present communication. 

Many unconditionally stable schemes have been proposed for the Cahn-Hilliard and related equations. The 
key ingredient of developing unconditionally stable schemes is appropriate temporal approximation of the chem¬ 
ical potential term, or the first derivative of the non-convex chemical free energy density with respect to chemical 
composition. Most often used are convex splitting methods 0E], where concave and convex parts of the chem¬ 
ical free energy are treated by different approximation methods. Convex splitting methods have been used in 
various related equations such as phase-field crystal equation [28], [25], the Cahn-Hilliard-Hele-Shaw system [27] . 
multicomponent Cahn-Hilliard equations [3 [23], and the Navier-Stokes-Cahn-Hilliard equation [10] [14]. The 
midpoint approximation method introduced in [6] is an alternative that does not introduce numerical dissipation 
and it has been used for the Cahn-Hilliard equation [5], the Allen-Cahn equation and the Cahn-Hilliard equation 
m, the Navier-Stokes-Cahn-Hilliard equation mm, liquid crystal dynamics m, and a two-phase flow model 
[T6] . The truncated Taylor expansion method was proposed for the Cahn-Hilliard equation in m The advan¬ 
tage of this method is that it can be extended to obtain unconditionally stable schemes for multicomponent 
Cahn-Hilliard equations CHI using truncated multivariate Taylor expansion. In these formulations remainder 
terms of the Taylor polynomials serve as numerical dissipation; indeed, if untruncated Taylor expansion were 
used for the standard single-variate non-convex quartic chemical free energy density, the resulting approximation 
would coincide with the midpoint approximation that has no numerical dissipation. Truncated Taylor expansion 
is also used in conjunction with convex splitting methods in 29, for the Cahn-Hilliard equation. Note that the 
composition gradient appearing in the Cahn-Hilliard equation is separately approximated in the above midpoint 
approximation method and Taylor expansion methods, but is treated as a unified contribution in convex splitting 
methods. Most of the above works consider quartic chemical free energy density. In m unconditionally stable 
schemes were developed for the Cahn-Hilliard equation for logarithmic free energy densities, using dedicated 
quadrature formula. 

The mechano-chemical free energy densities that we focus on in this work can be separated into two parts; 
viz. a logarithmic chemical free energy density used in [12] and a non-convex multivariate polynomial function 
that is eighth-order in strain and second-order in chemical composition, composition gradient, and strain gra¬ 
dient. Since the coupling of these terms makes the non-convex contribution rather complex, it is beneficial to 
develop a method that deals with components and component gradients in a unified framework instead of treat¬ 
ing each of them separately. Especially, a framework that provides accurate schemes in the presence of strain 
and strain gradients is crucial. To this end, we employ multivariate Taylor expansions as in [18] , but we now 
regard composition gradients and strain gradients as direct variables in addition to the composition and strain 
themselves, and use untruncated multivariate Taylor expansions to ensure unconditional stability and absence 
of numerical dissipation. Finally, the logarithmic chemical free energy density is treated in a similar way using 
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the Taylor expansion, partly drawing from the analysis in H2]. The resulting scheme also enjoys second-order 
accuracy as those proposed in [I2J \S \• To our knowledge this is the first unconditionally stable time-integration 
algorithm for mechano-chemical phase-transformation problems. Of note is the incorporation of gradient elas¬ 
ticity at finite strains—a feature that adds considerable complexity to the equations. The proposed method 
can be applied to any free energy density functions that are multivariate polynomials of components, compo¬ 
nent gradients, and higher-order component gradients to obtain unconditionally stable, second-order schemes. 
Finally, in addition to that they are accurate, these formulations allow for a straightforward simplification to 
produce practically useful second-order schemes that require less computation; those simplified formulations are 
called reduced formulations and are also investigated in this work. 

In Sec. [2] we present the variational formulation of our mechano-chemical problem. The corresponding 
spatially and temporally discrete formulations are developed in Sec. [3] Unconditional stability and second- 
order accuracy of our fully discrete formulation are studied in Sec. Finally, a suite of numerical examples is 
presented in Sec. [5] to demonstrate the performance of the algorithms in three dimensions. Final remarks are 
made and future work proposed in Sec. [6] 


2 Variational formulations for mechano-chemical spinodal de¬ 
composition 

We consider mechano-chemical spinodal decomposition in a body that occupies, at the initial time, a bounded 
domain Q in three-dimensional Euclidean space, in which we introduce the rectangular Cartesian coordinate 
system with Xj (J = 1,2,3) the corresponding coordinate variables. 

We are interested in the chemical composition field c(X,t) £ (0,1) and the mechanical displacement field 
u(X,t) in U. In this section we assume that these quantities and their spatial derivatives are continuously 
defined in U. The boundary of Q is assumed to be decomposed into a finite number of smooth surfaces T t , 
smooth curves Y t , and points so that dfl — TUTUS where T — U t r t , T — U t Y fc , and S =s Each 

surface V L and curve Y t is further divided into mutually exclusive Dirichlet and Neumann subsets that are 
represented, respectively, by superscripts of lowercase letters u , m, and g and those of uppercase letters T, 
M, and G, as I\ = r“ U if = If U rf and T t = T? U Tp. We also denote by T u = U t r“, r T = U t lf, 
r m = U,r"\ r M = U/Xf 7 , T 9 = U t Tf , and Y g = U,Yf ; the unions of the Dirichlet and Neumann boundaries. 
Our formulation and analysis presented in the following can be readily extended to include mixed boundary 
conditions. As in m, coordinate derivatives of a scalar function 0 are decomposed on T into normal and 
tangential components as: 


0, j = D(j)Nj + Djcj), 


where 


Dcj) := 4>,kN k , 

Dj<$> := (j),j - (P,kN k Nj, 

where Nj are the components of the unit outward normal to T. Here as elsewhere (-),j denotes the spatial 
derivative with respect to the reference coordinate variable Xj. 

Dirichlet boundary conditions for the displacement field u can now be given as: 

Ui — Ui on T u , Dm = fhi on T m , m — gi on T 9 , (1) 

where th, fhi, and gi are components of known vector functions on T u , T m , and T 9 . On the other hand, we 
denote the components of the standard surface traction on T t , the higher-order traction on T M , and the line 
traction on T G by ?*, Mi , and Gi , whose mathematical formulas will be clarified in Sec 
composition field c is assumed to have no Dirichlet boundary conditions throughout <9Q. 

2.1 Free energy 

We derive the initial and boundary value problem for mechano-chemical spinodal decomposition guided by 
variational consideration. The total free energy that we consider in this work is a functional of c and u defined 
as: 

n [c, u] := f T c + + T e dV- [ mTi d S- f D Ui Mi d S - [ Ui Gi dC, (2) 

Jn J r T J r M J r G 

where T c (c), T s (c, c a), and T e (c, Fu , Fu^k) are the chemical, surface, and mechanical free energy densities that 
are functions of chemical composition, c, gradient of chemical composition, c a, deformation gradient, Fu, and 


2.3 
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Figure 1: Plots of \Et (c G (0,1)) and contour plots of T e on the ~ e 3 plane at values of c corresponding to the 
maximum and the two minima of for a select set of parameters. Stable cubic/tetragonal crystal structures are 
also depicted. 


gradient of deformation gradient, Fu 7 k, at each fixed point XeO. Here, Fu — 8 u + Ui : j are the components 
of the deformation gradient tensor. These free-energy densities are explicitly defined as: 

T c := Hi (clog c + (1 - c) log (1 - c)) + H 2 c(l - c), ( 3 a ) 

:= ~c,aKab(c)c,b, ( 3 b ) 

T e := 81(c) (ei - e che m (c)) 2 

+ 82(c) (e 2 + e 2 ) + -83(0)03 (e 2 — 3 e 2 ) + 84(c) (+ e 2 ) + 85(c) (e 2 + e 2 + eg) 

+ 8 6 (c)(e 2,i + e 2 ;2 + e 2 ,3 + e 3,i + e 3,2 + e 3,3)> ( 3 c) 

where Hi and H2 are positive constants, e c hem(c), 81(c),..., 87(c) are polynomial functions of c, among which 81, 
84, 85, and Bq are positive, Kab(c) are also polynomial functions representing the components of the positive- 
definite, symmetric tensor that penalizes composition gradients, and finally ei,..., e6 are reparameterized strains 
defined as: 


ei = (8n + 822 + 8 33 )/>/ 3 , 

(4a) 

e2 = (8n — E 22 )/y/2 , 

(4b) 

e3 = (8n + 822 — 2 E 33 )/ V 6 , 

(4c) 

e4 — 823 = 832, 

(4d) 

65 — 813 — 831, 

(4e) 

e6 = 812 = 821, 

(4f) 


where Eu — 1/2 (FkiFkj — Sij) are the components of the Green-Lagrange strain tensor. Fig. [l] shows, for 
a select set of parameters, plots of T c (c G (0,1)) and contour plots of T e on the e 2 — e 3 plane at values 
of c corresponding to the maximum and the two minima of T c along with stable cubic/tetragonal crystal 
structures. To facilitate formulation and analysis, we denote by £ an array of c, c a, 8*j, and Fij^K and define 
a new multivariate scalar function T s + e := T s + T e . Note that, from definitions ( |3b| ), pc| , and Q, * S +e 
is a multivariate polynomial function of c, c a, Fu , and Fu^k, which is crucial for developing our accurate 
time-integration algorithm. 


2.2 Chemical equilibrium/non-equilibrium 

In this section we derive the equations for non-equilibrium chemistry. We first take the variational derivative of 
the total free-energy functional with respect to the chemical composition c in the direction of q to obtain: 


8 c U[c,u] = —n[c + eq, u] 


: f q(fl(c) + H(C))+q,AWA(C) dv, 
Jn 


(5) 
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where //(c), if(£), and Wa(C) are defined as: 


h 

dT c 

dc 

(6) 

H 

dc 

(7) 

W A 

d'I's+e 

dc,A 

(8) 


where //(c) is known as the homogeneous chemical potential. At equilibrium one has S c Il[c , it] = 0 and therefore 
from § one obtains the chemical-equilibrium equation as: 


/ 

Jn 


q (p(c) + H( 0) + q,AW A { C) dV = 0. 
We then apply the divergence theorem to Eqn. © to obtain: 


(9) 


/ 

Jn 


+ H- Wa,a ) dV + / q\V A N A dS = 0. 


L 


Standard variational arguments then lead us to the following strong form for chemical equilibrium and the 
corresponding boundary condition: 


// + H — Wa,a = 0 in Q, 
WaNa = 0 on f. 


(10a) 

(10b) 


Note that one can further write Eqns. (10a) and (|10b| as: 


// + H — (Kabc,b),a —— 0 in Q, 

(Kabc j b)Na 0 on r, 

but in this work we avoid this step to simplify our stability analysis presented in Sec |4.2| We identify the 
left-hand side of Eqn. (10a) as the chemical potential //, i.e.: 

// = // + H - Wa,a • 


(ii) 


With the expression for the chemical potential (11) in hand, one can formulate the non-equilibrium chemistry 
problem using the mass balance law in conjunction with the phenomenological representation of the flux, that 
is: 


Dc . 

Wt +jA,A-°, 


( 12 ) 


where, as elsewhere, D/D/ represents the material time-derivative, and the flux is in coordinate notation, 

Ja — —LabH,b, (13) 

with Lab(c) being the components of the positive definite mobility tensor. Eqn.(|12|) requires another boundary 
condition for JaNa on T. Throughout this work we follow the proposition in [12) and set: 

JaN a = 0 on T. (14) 

In this work we adopt the mixed formulation in which // as well as c are regarded as primary unknowns. 
Multiplying Eqns. (12) and (11) by admissible test functions q and zq applying the divergence theorem applying 
boundary conditions (10b) and (14), one obtains the weak form for the non-equilibrium chemistry problem as: 

(15a) 


/.(■ 
L 


q ]] • + q,AL AB {c)n,B ) dV = 0, 


(«/ (-A* + P(c) + H{ 0) + v, a Wa{ 0) dv = 0. 


(15b) 
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2.3 Mechanical equilibrium 

To formulate the problem of mechanical equilibrium, we take the variational derivative of the total free energy 
§ with respect to u that satisfies the Dirichlet boundary conditions ([!]). The test function w is then to satisfy: 

Wi — 0 on T n , Dwi = 0 on T m , Wi = 0 on T 9 . 

The variational derivative with respect to u is then obtained as: 


S u U\c. u] — ——II[c, u + ew] 
de 


= [ (Wi,jPij(C) + Wi^JKBiJK ( 0 ) d V - f WiTi d S - [ DwiMi d S - [ WiGi AC, 
J n Jr T Jr M Jt g 


(16) 


where Pij(C) are the components of the first Piola-Kirchhoff stress tensor and Buk^C) are the components of 
the higher-order stress tensor that are defined as: 

d^s+e 


Pu := 

BiJK ■ = 


dFu ’ 

d^s+e 


dFi 


iJ,K 


At equilibrium one has 5 u Tl[c, u] = 0 so that from (16) one obtains: 

[ (wi,jPij(C) + Wi,jKBijK( 0) d V - [ WiTi d S - [ DwiMi d S- f Wi Gi d C = 0. (17) 

Jn Jv T Jr M J r G 


Eqn. (17) is the weak form for mechanical equilibrium that is to be solved in conjunction with the weak form 


for the non-equilibrium problem of chemistry (15) 


The variational argument can further lead us to identify the strong form and the Neumann boundary con- 

in Q. 


ditions corresponding to (17) as the following: 

PiJ,J — BiJK,JK — 0 

PijNj — BijK,i<Nj — Dj^BijKNx) + BijK ( BllNjNk — bjx) — T 

BijKNxNj — Mi 

lBi JK N K N V j\ = Gi 

where bu are the components of the second fundamental form on T t , Nj are the components of the unit outward 
normal to the boundary curve T t C T t /, and, on each Yf, {BukNkNj ] := B iJK N+N t + + BukNkN*- is 
the jump, where superscripts + and — represent two surfaces sharing Yp; see m for details. 


on T 1 , 

p AT 

on 1 , 

'Y' G 

on Y , 


3 Numerical formulations 

3.1 Spatial discretization 


We now discretize Eqns. (15) and (17) in space for formulations that are amenable to numerical analysis. 


Since T s+e is a function of Pij,K , our weak forms (15) and (17) involve second-order spatial derivatives of the 
displacement field u and the test function w : and thus require them to be W 2 ’ 2 , where IV s,p is the standard 
Sobolev space. We denote by S h an appropriate finite-dimensional subspace of W 2,2 (fl) and define: 

Vu — ^u h G [S h ] 3 : Ui — Ui on T w , Dul — fhi on T m , iq = <ji on Y^| , 

G [S h ] 3 : = 0 on T n , Dw^ = 0 on T m , w*t = 0 on Y^| , 

assuming that S h allows for exact representation of the Dirichlet boundary conditions ([!]). We also denote by 
T h an appropriate finite-dimensional subspace of W 1,2 (fl). 

The space-discrete counterparts of the weak formulations ( |15| and © are then given as the following: 

Seek c h (X,t) E T h x [0 ,T\, p h {X,t) E T h x [0,T], u h (X,t) E Vp x [0 ,T] such that for all g h (X) E T h , 

v h {X) E T h , w h (X) E Vp: 

( gh ^f + A WAA) dV = 0, 

J' P (~H h + n{c h ) + H(C h )) + v*.AW A (C h )) dV = 0, 

I Wi,jPij{C. h ) + w^ jk Buk( c' 1 )) AV - J wifi AS - j 


J 

Jn 


(18a) 

(18b) 


p M 


DwiMi AS- f wlGi AC = 0, (18c) 

J r G 
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where £ h (X,t) is an array of c h , c h Al F^j(— 5ij + u^j), and F^ JK . Here, all spatial derivatives are now to be 
understood in the weak senseQ 

Note that, setting q h = 1, v h = 0, and w h = 0 in the weak form (18), one recovers the mass conservation 
law: 

/ c 1 dV = 0. 


d t 


L 


On the other hand, assuming that all Dirichlet and Neumann boundary conditions are time-independent and 
setting q h m v h — Dc h /Dt, and w h = Du h /Dt in Eqns. (18) and adding them together, one obtains: 


_drr 

dt 


- [ ^aLab^b dV, 

Jn 


(19) 


where U h is the space-discrete total free energy at arbitrary time t defined as: 

U h (t)= [ f^ c (c h ) + ^ a+ e(C h )) dV - [ u^TidS- [ Du^MidS- f u!?Gi dC, ( 20 ) 

Jn v J Jv T Jr M Jr G 

spatial derivatives being understood in the weak sense. Recalling that the mobility tensor is positive definite, 
Eqn.(19) implies non-increasing free energy. 

These two properties, mass conservation and non-increasing total free energy, are to be inherited by our 
space-time discrete formulation developed in Sec. |3.2[ specifically, the latter property furnishes the notion of 
stability. 

As is well known, substituting Eqns. (13), (11) and ([6j{8]) into (12) leads to a fourth-order PDE in strong 
form. Its weak counterpart has up to second-order spatial derivatives on the composition, and was the basis 
for Discontinuous Galerkin-based finite element methods in the work of Wells et al. [23 > a nd m ore recently 
for C 1 -continuous IGA-based methods in EH- Here, we use the split formulation of Eqn. ([15]), because, as 
shown above, it permits the fields c h and fi h to be chosen to lie in the same space, T h , in the resulting finite 
dimensional statement of the full problem ( |18| . This coincidence of spaces is crucial for satisfaction of the 
fundamental stability result just derived, and for its extension to the time-discrete setting. 

3.2 Temporal discretization 


We proceed to discretize Eqns. (18) in time to obtain a formulation that produces a solution at time t n+1 given 
a solution at time t n . The proposed time-discrete formulation is given as the following: 

Given c h ’ n (X) G T\ //’ n (X ) e T\ u h ’ n (X) e Vi, seek c h ’ n+1 (X) e T\ /a h ’ n+1 (X) e T h , u h ' n+1 (X) e Vi 


such that for all q h (X) e T , ^(X) G T , w h (X) e V^: 
^ n 


(21a) 

(21b) 


i if {%} + dV = 0, 

J {y h (-{mT + {n{c h )Y + (ff(C' l )} n ) + v h A{w A {C, h )} n ) dV = o, 

J (™?,j{Pu(C h )} n +wZjK{BuK(<: h )} n ) dV 

-[ wUfi} n dS-[ Dw? {Mi\ n dS - [ Wi {Gi] n dC = 0, (21c) 

J r T J r M Jt g 

where terms with braces {#} n (X) represent time-discretizations of those quantities inside the braces • (X,t) 
on the time-interval t G [t n , t n+1 ], and they are defined in the rest of this section so that the formulation in ( |21| 
is second-order accurate and unconditionally stable. 

The stability analysis presented in Sec. 


is greatly motivated by that presented in T2j for the Cahn- 
Hilliard equation, where {j2(c h )} n was defined using dedicated quadrature formulas so that a non-increasing 
chemical free energy would be the direct consequence of the weak formulation. In this work we follow the same 
course, but, instead of developing special quadrature formulas, we employ Taylor expansions, which enables one 
to simplify the argument as well as to extend the stability analysis to coupling with gradient elasticity. 

For convenience we denote by C h,n (X) the temporal approximation to £^(X,t n ). In addition we define 

Af — y.n+1 _.n \ h — h,n+l_h,n A h — h,n+l _ h,n a j?h — h,n+l _ h,n i a T? h — 

Lai .— t t , Zac •— c c , zac ^4 .— c A c A , ZA r iJ .— u i J u i J , ana LAr iJ K .— u i JK u i,jk- 


1 Since S h C yV 2,p (H), second spatial derivatives of w h and u h are properly defined only in a weak sense—a technicality that is 
usually not emphasized in finite dimensional weak forms. 
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We first represent {/ i(c h )} n in terms of c h,n and c h,n+1 at each fixed point X £ Q. To this end we observe 
that the Taylor expansion of T c (c) around c h,n+1 leads to the following identity: 


y c (c h ’ n ) = V c (c h ’ n+1 ) - jl(c h ’ n+1 )Ac h + \^(c h ’ n+1 ){Ac h f - \^{c h ’ n+1 ){Ac h ) 3 + 

2 dc o ac z 24 dc d 


= V c (c h ’ n+1 ) - (W’ n+1 ) - \^{c h ’ n+1 )Ac h + iTf(c‘ 
V 2 dc o ac z 


h , n -\-1 


)(Ac 


Ac' 1 + ^ 7 4^(e(Ac' 1 ) 4 , 


24 dc 3 


( 22 ) 


where £ = (1 — a)c h,n + ac h,n+1 for some a (0 < a < 1) by Taylor’s Remainder Theorem. We then define 
{fl(c h )} n as the quantity in the parentheses in |22|, that is: 


{fi(c h )} n := jl(c h ’ n+1 ) - ±2^(c h ’ n+1 )Ac h + ^(c h ’ n+1 

Z CLG O QC 


)(A c h )\ 


(23) 


so that: 


{fi(c h )} n Ac h = V c (c h ' n+1 ) - '5> c (c h ’ n ) + 


1 d 3 /x 
24 “d^ 


(0 (A c h 


(24) 


at each fixed point X £ Q. Identity ( |24| becomes a convenient tool in the stability analysis encountered in Sec. 
[ 4 ] noting especially that d 3 /i/dc 3 (£) > 0 by virtue of ( |3a| . 

We dehne {if(C^)} n , {WA(C h )} n 1 {F > ij(C /l )} n , and { BijK(C h )} n in a similar fashion. We first denote by 
V [0; ft c , ftvc, ff, fvf] the function obtained by applying operators (d/dc)Ac h , (<9/<9c a) AcJa, (d/dFij)AFj}r, 
and K respectively ft c , ^v c , ^f, and fvf (ft c , fvc, ^f, fvf > 0) times to a scalar-valued 

multivariate function 0(C)- For instance we have: 


£>[0;O,O,O,O] = 0, 

^[0; 1,0, 2,1] = 


0 4 0" 


dcdFijdFkLdFmN,o 


A^AF^jAF^AF^o 


= V 


d<p 

dF~j 


; 1 , 0 , 1,1 


AF?j. 


We also define f = k c + fvc + ^f + fvf- The Taylor expansion of T s+e (C) around £ h,n (X) at a fixed point 
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X G O then leads to the following identity: 


*s±e(C h ’ n+1 ) ^s+e(C h ’ n ) + J2 

K>1 

= Vs + e(C h ' n ) + E 


Av c > 1 
/“C> 1 


+ E 


ft! 

Fc_1_ 

k k c !«vc^f!«vf! 

ftVc 1 


P [T s+e ; «c, ^Vc, «F, ^vf] (C^’ n ) 


P [T s+e ; «c, ftVc, «F, ^Vf] (C^’™) 


«Vc>! 

K> 1 


ft ft c !ft Vc*^F F • 


P [T s+e ; FC c , ftVc, «F, ^Vf] 


+ E 

kj?M 

Av> 1 

+ E 


fti? 


! Kl ,.! K „l to / [$s+e; Kc ’ KVc ’ Kf ’ Kvf ' (C ’ n) 

P [^ s + e ; «c, Kvc, kf, kvf] (^ ,l ’ n ) 


ft ft c !ftVc^F'^VF' 
ft V F 1 


^ V F >1 

K> 1 


ft ft c !ft Vc F !ftVF • 


/ 


= ^ S +e(C' l ’ n ) + 


+ 


+ 


+ 


E -- 

E^ K- ( 


' ft (ftc — 1)!ftVc'Fir 1 !ftVF! 

I 

V 

t 


E A - 

E—/ /C / 


ft ft c ! (ftVc — ldftidftVF! 

\ KVc>! 

\ «>1 


V [H ; ft c — 1, ftvc, ff, fvf] (C^’ n ) I ^° h 


V [Wa; k c , kvc — 1, ff, fvf] (C^ ,n ) Ac^a 


/ 


E 1 - 

E—/ tc / 


1 


\ 


^ ft ft c !ftv c ! (ff - 1)!fvf! 

y K F >1 v 7 

\ K>1 


P [Pjj; K c , KVc KF — 1 , kvf] (£ h ’ n ) 


A Fj 


/ 


/ 


E A - 

E—/ K / 


\ 


_ ft ft c ! k Vc IkfI (kvf - 1)! 

I k V f> 1 

\ K>1 


[BiJic; ftc, ftvc, ftF, FVF — 1] (C^ n ) 


) 


AF iJK , 

(25) 


where summations are over all possible combinations of ft c , ftvc, ftF, and ftvF for each ft. These summations 
are finite as T s + e (C) is a multivariate polynomial function of c, c a, Fu, and Fu,k- Factors in the summation 
in the first line arise since V [T s+e ; ft c , ^v c , ff, ^vf] ( C h,n ) appears in a straightforward Taylor-series expansion 
ft!/ft c !ftv c !FF!FvF! fi mes ^ue to this number of possible permutations; for the sufficiently smooth , T considered 
here, the following terms all reduce to (1/3!) • V [T s + e ; 0, 0, 2,1] (C^ ,n ) and therefore this term in the above 
summation is to be multiplied by 3!/0!0!2!l! = 3: 


1 d 3 V s+e 

3! dF u dF kL dF mN ,o 
1 d 3 V a+e 

3! dF u dFmN,odF kL 

1 d 3 V a+e 

3! dF mN ,odF u dF kL 


(C; h ' rl )AFijAF kL AF rnNt o, 

(C h ' n )AF u AF mN , 0 AF kL , 

(C h ' n )AF mN ,oAFuAF kL . 


We then define {H(£ h )} n , {WA(C, h )} n , {P,j(C,' )}”, and { B,.j k(C h )}" as those quantities in the parentheses in 
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([25]), or: 

{H(C h )} n 

■■= H(C h n + l (^-(C h ’ n )Ac h + |^(C' l ’ n ) Ac% + ^L(C h ’ n )A ff M + ^—(C^AF^n) + R c (C h ' n ), 

I \ OC OC,B OriM OriM,N J 

(26a) 

{W A {C h )} n 

TJ7 / >Ai,ro\ , 1 / 9Wa / j,h ,n\ a h , &WA/±h,n\ \ h , a nh , 9Wa / >/i,n \ a t-i/i \ , D Vc/>(i,n\ 

:=Wa(c ) + 2 ^ (C } + ^y (c )Ac ’ s + ^w (c )AF(M + aw (c )A H + i(C ); 


(26b) 


{P,J« h )} n 

:= Pi J (C h,n ) + ^ (j^(C h ’ n )&c h + ^ " h 


{iWK' 1 )}" 


(C’ n )Ac» + ^(C’ n )Ac% + ^-(C h ’ n )AF ? m + i^+C' 71 ) A^Vat ) + tfjK'*'”), 


0P,. 


dFiM,? 


(26c) 


: = 


where 


/An , 1 fdB iJK / M,h,n\ a h , a# 

iJK ,+h,n\ a h dB iJK ,j.h,n\ a 77 A 1 iJK ,^h,n\ * rih \ . D VF /+h,n\ 

:(C ’ ) + « ~ (C )Ac + -(C )Ac )B + —^— (c )AF ZM + —-(C )AF ZMjiV + iW(C ), 

^ V °^ c s OtiM or im,N J 

(26d) 


i? c : » V - t -—j-j-j--P [ 17 ; — 1 , kvc, «vf] 

« (« c ■ — IjIkvJkfIkvf! 

K c A_ 1 


flX C := £ 1 

^' K 


K\7c^ 

k>3 


1 k kJ (^Vc - 


22 [VFa; ac c , acvc — 1 , acf, a^vf] ; 


: = E 1 

Z —' K, 


KJ7 > 

k>3 


i ft ft c !ftvJ (ftF — 1)!acvf! 


nVF . _ 1 

n i.jK • — y — 

' K 


F ^ 
k>3 


i ft k c !kvc^f! (ftvF — 1)! 


P [Pij; ft c , A^Vc, A^F — 1, A€vf] ; 


22 ft c , ftVc, a^F, ftVF — 1] ; 


(27a) 

(27b) 

(27c) 

(27d) 


so that: 


at each fixed point X £ Q. Finally, other quantities in Eqns.(21) are defined as: 

(28) 

f Dc h \ n 

A c h 

(29a) 

{-Dt} : 

~At’ 

{c h } n 

= [c h ' n+1 + c h ' n ]/2, 

(29b) 

{p h r 

= [,u h ’ n+1 + n h ' n }/2, 

(29c) 

{LAB(c h )} n 

= [L AB (c h ’ n+1 ) + L AB (c h ’ n )]/ 2, 

(29d) 

m n 

= [T" +1 + Tf ]/2, 

(29e) 

{Mi} n 

= [Mf +1 + M?]/ 2, 

(29f) 

m n 

= [G( ,+1 +Gr]/2, 

(29g) 

where T™(X), M z n (X), and G™(X) are the components of the boundary tractions at t n , Ti(X,t n ] 

), Mi(X,t n ), 


and G;(X,t n ), respectively. 


4 Analysis 

In this section we prove mass conservation, unconditional stability, and second-order accuracy of the time- 
integration algorithm proposed in Sec |3.2| 
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4.1 Mass conservation 


Provided that Eqns.(21| are satisfied for all q h (X) 6 T h , v h {X) € T h , and w h (X) e V^,, those equations are 
necessarily satisfied when we substitute the following for these test functions: 

q = 1, v h = 0, w h = 0. 


One then readily obtains: 

[ c h ’ n+1 dV= f c ' n dV, 
Jn Jn 

which implies that mass is conserved from time t n to time t n+1 . 


4.2 Stability 

We now investigate stability of the proposed time-integration algorithm. To this end, we assume that all Dirichlet 
and Neumann boundary conditions are time-independent; that is Ui, mi, gi, Ti , Mi, and Gi are constant in 
time. 

Provided that Eqns. (21) are satisfied for all q h (X) G T h , v h (X) G T h , and w h (X) G V^;, those equations 
are necessarily satisfied when we substitute the following for these test functions: 


h r h^i n 

q = I , 


h 

V — 


_ /yb\ 

\Dt] ’ 


h 

w = 


At 


We then add the resulting three equations together and use identities (24) and ( |28| to obtain: 


At 


■f a (i» h }y{LAB(c h )r{» h }: B + dy ’ (30) 


where II^ ,n is the space-time discrete total free energy at t = t n defined as: 


it 


: [ T c (c /l ’ n ) + T s+e (C /l ’ ri ) dE- 

Jn 


f u^ n Ti d S - [ Du^ n Mi d S - 

JrT Jr 


J pM 


/ 


dC, 


(31) 


where spatial derivatives should be understood in the weak sense. Note that, as the mobility tensor is positive 
definite and d 3 /i/dc 3 (£) is positive by definition (3a), the right-hand side of (30) is non-positive. Eqn. (30) 
therefore states that the discrete total free energy is non-increasing from time t n to time t n+1 and the algorithm 


proposed in Sec |3.2| is necessarily unconditionally stable. It should also be noted that, as seen in Eqn. (30), 
numerical dissipation originates only from approximation of the logarithmic chemical potential {g(c h )} n and not 
from approximations {H(£ h )} n , an d {B iJK(C h )} n 5 our multivariate Taylor expansion 

method applied to a multivariate polynomial free energy produces no numerical dissipation. 


4.3 Consistency and second-order accuracy 

We proceed to show second-order accuracy of the proposed scheme. Following the standard treatment for the 
consistency analysis, we replace c h,n (X ), fi h,n (X ), u h,n (X), c h,n+1 (X), fi h,n+1 (X), and u h,n+1 (X) in the time- 
discrete formulations (21a), (21b), and (21c) with the corresponding solutions to the time-continuous problem 
(18) at t n and t n+1 ; we denote the left-hand sides of the resulting equations by 7™, 7™, and 7™, respectively. The 
following approximations for an arbitrary function <p(X,t) readily obtained by Taylor expansions are of use: 

«xx+')<-Hx, n = , {X<n+ Am (x , n+ 


HX,f 


,n+l\ 


■4>(x,t” 


At 


= ^(x,n++o(At 2 ), 

Dcj> 


A0(X) ■.= ^(X,t n+1 )-4>(X,t n )=^(X,t n )At + 0(At 2 


The definitions of the high-order terms R c , R ^ c , Rfj , and RJj F k in (27) show them to be 0(At 2 ). Therefore, one 
can, for instance, show that {Pij(C h )} n given in (26c) can after substitution of the time-continuous solutions 
be approximated as: 


{Pu(C h )} n = Pij(X, t n ) + + 0(At 2 ), 
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where Pu(X,t) := Pij (^ h (X,t)). Treating other terms similarly, we can readily show the following: 


ic = ic(t n ) + nn+o(At 2 ), 

i; = iAt n ) + ^^(t n ) + o(At 2 ), 

In = Iu{t n ) + + 0{At2) ’ 


(32a) 

(32b) 

(32c) 


where / c (£), (£), and I u (t) are, respectively, the left-hand sides of Eqns. (18a), (18b), and (18c). Since 


the first two terms of the right-hand side of each equation in (32) are zero, one concludes that the proposed 
time-integration scheme ( |21| is of order 2. 

We note here that in the above consistency analysis the specific formulas for R c : Ra°i and R^j F k are 
unimportant. Indeed one can ignore some or all high-order terms existing in R c , Ra°, Rfj , and RYjk when 


evaluating (26), with the resulting reduced formulations remaining second-order accurate. Such reduced formu¬ 
lations lose unconditional stability, but are often equipped with practical accuracy. Requiring less computation, 
they can serve as good alternatives to the original formulation in many problems; see Sec|5.3|for examples. 


5 Numerical examples 

In this section we demonstrate the robustness and accuracy of the numerical formulation proposed in Sec. [3] 
Parameters for the chemical free energy density function T c introduced in Eqn. (3a) are set as A\ — 1 and 


^4.2 — 3 so that this function takes its maximum at c— 0.5000 and minimum around c = 0.0707 and 0.9293. 

Fig. [2] shows plots of T c against chemical composition c, where the chemical spinodal region is also shown 
as the interval between two dashed lines. We denote by d{— 0.1152) the difference between the maximum and 
minimum values of T c . The diffusion and mobility tensors appearing in Eqn. ( |3b| ) and in Eqn. ( |12| ), respectively, 
are defined as Kab(c) — 2 ~ 6 5ab and Lab(c) = 6c(l — c)5ab- Parameters for the mechanical free energy density 
function T e defined in Eqn. pel are given as e c hem = — l/16(c —0.48), £>i(c) — 13/4r, £> 2 ( 0 ) = —5/32r(c —0.50), 


£?3(c) = l/4r(c — 0.50), B^c) = r , £>s(c) = 2r, and Bq(c ) =2 -14 rc, where r = 2 11 d. These free energy density 
function parameters have been chosen to produce a suitable microstructure; a wide range of microstructures 
results by variation of these parameters. The free energy density component T e thus defined characterizes 
crystallographic structural changes between a cubic phase and three energetically-equivalent tetragonal phases; 
the cubic phase loses stability and transforms into one of three stable tetragonal phases as the local composition 
c increases. This transformation is represented by the projection of T e onto the e 2 — e 3 plane, showing a 
continuous transition from convex to non-convex form with three-wells on the e 2 — e 3 plane for c > 0.5; see Fig. 
[3] The distance on the e 2 — e 3 plane from the origin to these three minima is designed to be 1/4 when c = 1; 
that is, minima occur at (+\/3/8,+1/8), (—\/3/8, +1/8), and (0, —1/4), each corresponding to a tetragonal 
crystal structure that is elongated in the Xi-, X2-, or X3- direction, respectively. Stable crystal structures are 
also depicted in Fig. [3] 

We also note here that, to compute {H(£ h )} n , {WA(C h )} U i {Pij(C h )} n i an d {BijK(C h )} n defined in (26) 
and (|27|), one needs to take derivatives of T s + e (C) with respect to c, ca, P%j , and Fij,k , respectively, 2, 2, 8, and 



Figure 2: The chemical free energy density function T c , plotted against chemical composition c for parameters 
used in our example problems. 
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(a) (b) (c) 

Figure 3: Contour plots of the projection of the mechanical free energy density function 4/ e onto the e 2 — e% plane 
for (Ja| c = 0.0707, Q c = 0.5000, and ([c]) c = 0.9293. The corresponding stable crystal structures are also depicted. 


2 times. This follows from Equations (3a)— ([3c]) and (|4a|- It is therefore sufficient to set k c < 2, kvc < 2, 

< 8, and kvf < 2. 

The reference domain occupies a unit cube Q = (0, l) 3 and we set boundary conditions as displacement 
Ui = 0 and displacement normal gradient Dm — 0 on X3 = 0, traction Ti — 0 and higher-order traction Mi — 0 
on X3 = 1, and normal displacement UkNk = 0, shear traction Ti — T^N^Ni — 0, and higher-order traction 
Mi = 0 on Xi,X 2 = {0,1}. Further, we set the higher-order traction jump Gi = 0 on all edges where m = 0 is 
not specified. 

In this work we employ isogeometric analysis (IGA) that allows for the easy construction of C p -continuous 
basis functions for arbitrary order p. IGA has previously been adopted to treat higher-order spatial derivatives, 
e.g., in 11, 20 j and we refer the reader to these works for details on the development of IGA basis functions 
with the desired degree of continuity. Our three-dimensional IGA basis functions are formed by the tensor 
products of one-dimensional, second-order B-spline basis functions on uniformly spaced knots. The initial 
condition of the local composition c is produced on a coarse 2 3 mesh that has 2 elements , or 4 basis functions, 
in each direction. The primitive one-dimensional B-spline basis functions in the X±-, X2- , and X 3 -directions 
are indexed as i± ,Z 2 , and 23 , where AG 2 G 3 = 0,1,2,3, and control points for A, 22,23 = 1,2 are given as 
0.48 + 0.01 sin(999sin(997ii + 99122 + 98323 + 1)) and those for 21 , 22,23 = 0,3 are computed to satisfy the 
boundary condition (10b). This initial condition is projected onto finer meshes (16 3 , 32 3 , and 64 3 meshes) 
exactly by successive uniform h-refinements by knot-insertions [3]. Finally, initial conditions for the chemical 
potential and the displacement field are given as /x = 0 and m = 0. 

We solve these three-dimensional, mechano-chemical, initial and boundary value problems using the uncondi¬ 
tionally stable, second-order accurate time-integration algorithm proposed in Sec[3]with absolute residual error 
tolerance set to 10 -10 . Our open source research code [22] is written in C, uses PETSc 3.5 for linear/nonlinear 
solvers and mathgl 2.3 for plotting. We also utilized Mathematica 10 to compute high-order indicial sums 
appearing in (26) and (27), and tangent matrices required for nonlinear solvers. 


5.1 Temporal accuracy 

In this section we study the temporal accuracy of our proposed formulation using a 32 s mesh throughout. We 
use progressively finer timesteps At = 4x 10 -3 , 2x 10 -3 , 1 x 10 -3 , 5x 10 —4 , and 2.5 x 10 -4 and integrate up to 
t = 4, at which time the solutions were found to be almost at steady state. 

Fig. [JJshows color plots of e 2 along with contour curves of c for solutions obtained using At = 4xl0 -3 , 2xl0 -3 , 
1 x 10 -3 , and 5x 10 -4 . While At = 4x 10 -3 leads to a completely different morphological evolution implying 
insufficient temporal resolution, the absence of visible differences between the solutions for At = 2xl0 -3 , lxlO -3 , 
and 5x 10 -4 indicates convergence of the microstructure as timesteps are refined. 

Figs. [5j|a| and [HJJb]) show plots of discrete total free energy H h ’ n against time for t G [0,4] and t G [0.9,1.1], 
respectively, for solutions corresponding to At = 4xl0 -3 , 2xl0 -3 , lxlO -3 , and 5xl0 -4 . One can observe 
convergence of the free energies of the solutions on t G [0,4] with time step refinement; we also remark that, 
while the numerical solution with At = 4 X 10 -3 has a markedly different microstructure from those computed 
with finer timesteps (Fig. it does indeed converge to the same free energy at the steady state (Fig. [5[a])). 
Of note is that the discrete free energy is non-increasing for any time step, implying unconditional stability of 
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£ = 0.08 


£ = 0.10 


£ = 1.00 


£ = 2.00 


£ = 4.00 


A£ = 4 xlO -3 


A£ = 2 xlO -3 


A£ = 1 xlO -3 


A£ = 5 xlO -4 












Figure 4: Plots of e 2 and contour curves of c on deformed configurations at select times for solutions corresponding 
to four different time steps A£ = 4 x 10 -3 , 2 x 10 -3 , 1 x 10 -3 , and 5 x 10 -4 on a 32 s mesh. 
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(a) (b) (c) 

Figure 5: Plots of discrete total free energy U h,n against time t for ([a| t £ [0,4] and ([b]) t E [0.9,1.1] for four 
different timesteps At = 4xl0 -3 , 2xl0 -3 , lxlO -3 , and 5xl0 -4 . 0 Plots of error, \\c h,n — c h \\ 2 , at t = 1 against 
timestep At on a logarithmic scale. A fixed 32 3 mesh was used. 


our proposed time-discrete formulation as expected from the analysis in Sec. |4.2| In this regard we point out 
that, though the scheme is unconditionally stable, the time-step size still sets a constraint on the solvability of 
the nonlinear system; the nonlinear solver quit converging to the set tolerance around t = 0.088 when we used 
At = 8 x 10 3 , which is in any case too large to have a reasonable solution for this specific problem. 

Finally, we regard the solution for At = 2.5xl0 -4 as exact and compute L 2 -norm of the error of the solution 
field c, || c h ' n — c h || 2 , at t = 1 for each solution. Fig. jHJc) shows plots of error \\c h,n — c h || 9 against time step At 
on a logarithmic scale, where one observes a second-order temporal convergence as expected from the analysis 
in Sec. 14.31 

5.2 Spatial convergence 

We continue in this section to investigate spatial convergence of the solutions with mesh refinement. We solve 
problems on three different meshes, 16 3 , 32 s , and 64 3 , using a fixed timestep At — 2xl0 -3 that is regarded 
as small enough. Time-integration is again performed up to t = 4. Fig. [ 6 ] shows the temporal evolution of the 
microstructure on these three meshes. One observes that, while the 16 3 mesh seems under-resolved, no further 
morphological changes appear under refinement from 32 s mesh to 64 3 mesh, giving good evidence of spatial 
convergence of the microstructure. The solution on a 64 3 mesh with At = 1 X 10 -3 is also plotted to show that 
At = 2 xl 0 -3 gives sufficient temporal resolution for this spatial convergence analysis. Figs. [7J[a]) and[7J[b| show 
plots of corresponding discrete total free energy U. h,n against time for t £ [0,4] and t £ [0.9,1.1], respectively, 
where one can further observe convergence with respect to free energy of the numerical solutions under mesh 
refinement. Discrete total free energy plots for the solution on the 64 3 mesh with At — lxlO -3 are also shown, 
which provides further assurance that the timestep of At — 2 x 10 -3 is small enough for this observation. 

We conclude this section by investigating the converged microstructure obtained in our numerical analysis. 
Figs. | 8 j|a| and [ 8 j[b]) show the top views of color plots of and e 3 in the deformed configuration at t — 4 
computed on a 64x mesh with At = 2 x 10 -3 . For the sake of visualization of the deformation, these top 
views are overlaid with distorted 32 2 meshes. The deformation reveals twin boundaries between two of the 
three tetragonal variants, viz. those two living in the upper-half plane in Fig. [3j|c|. Figs. [ 8 j|a| and| 8 j|b| also 
make clear the large deformations that the microstructures suffer; note the distorted mesh corroborating the 
strain values in the legends of Figs. [4] and [ 6 ] Fig. | 8 j|c| shows dot plots of (e 2 ,es) for a square array of 64 x 64 
points, uniformly spaced along X\ and X 2 , and lying on the X3 — 1 plane. Also reproduced, are the contour 
plots originally shown in Fig. [3j|c]). The sharp localization of the strain state in e 2 — e% space is understood 
to be a consequence of the specific boundary conditions employed for this problem. We draw attention to the 
localization in the vicinity of only two of the three wells in e 2 — e% space. The points lying between the two 
wells are found to lie in the twin boundaries between the corresponding variants on the physical domain Q. In 
other computations we have found the equidistribution of the strain state across all three wells when perfectly 
symmetric boundary conditions are used. 
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£ = 0.10 
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Figure 6: Plots of e 2 and contour curves of c on deformed configurations at select times for solutions computed on 
three different meshes, 16 3 , 32 s and 64 3 with At = 2xl0 -3 . The solution on a 64 3 mesh with At = lxlO -3 is also 
shown for comparison. 




(a) (b) 

Figure 7: Plots of discrete total free energy H h ’ n against time £ for ([a| £ G [0,4] and (|b| £ G [0.9,1.1] for three 
different meshes, 16 3 , 32 s , and 64 3 , with timestep At = 2x 10 -3 . The solution on a 64 3 mesh with At = 1 x 10 -3 
is also shown for comparison. 
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(a) 


(b) 


(c) 


Figure 8: Plots of (a| e 2 and ([b]) on Z = 1 obtained on the 64 3 mesh with At = 2 x 10 —3 , overlaid with 32 2 
plotting meshes, dc) Dot plots of e 2 and e% values on X 3 = 1 evaluated over a uniformly spaced square array of 
64 x 64 points, laid over the contour plots originally shown in Fig. K3|jc|). 


5.3 Reduced formulations 


In the preceding numerical examples, we have evaluated {/i(c h )} n , {H(£ h )} n , {Wa(C h )} n > {Pij{C h )} n > and 
{B iJK(C h )} n following ( [26] ) while using the full summations in ( [27| given by k c < 2, k^/ c < 2, kf < 8, and 
kvf < 2. In this section we introduce reduced formulations in which some high-order terms in {fl(c h )} n , 
{H(£ h )} n , {WA(C h )} n ? {Pij(C h )} n ? an d {5ijK(C^)} n are ignored; specifically, we consider two formulations 
obtained by setting kf < 2 and < 4 instead of < 8. These reduced formulations are not unconditionally 
stable, but often provide solutions of sufficient quality at lower computational cost. 

To demonstrate this point we solve the example problem encountered in Sec. 5T on a 32 3 mesh with 
At = 5 x 10 4 up to t = 4 using these reduced formulations and compare them with the full formulation. Using 
2.60GHz Intel Xeon E5-2670 processors on 8 x 8 x 8 partitions, the actual time required for the time-integration 
was measured for each formulation. The results are summarized in Table [T] The reduced formulation with 
kf < 2 was about 2x as fast as the full formulation with kf < 8, but the solution diverged around t— 1, while 
the reduced formulation with kf < 4 was about 1.5x as fast, retaining practical stability. 

Fig. shows color plots of e 2 and contour curves of c for these solutions at select times t, where one can 
observe that microstructures obtained by the reduced formulations are almost identical to those obtained by 
the full formulation, although the reduced formulation with kf < 2 diverged around t — 1. 

Figs. [ 10 R and [lQfb| show plots of corresponding discrete total free energy H h ’ n against time for t £ [0,4] 
and t £ [0.9,1.1], respectively. The energy corresponding to the reduced formulation with kf < 2 takes slightly 
larger values than the other two formulations until it diverges in the neighborhood of t = 1. In contrast, 
the energy corresponding to the reduced formulation with kf < 4 is indistinguishable from that of the full 
formulation. 

This study indicates that, despite the possible loss of unconditional stability, reduced formulations can 
produce solutions that are of sufficient quality in practice and can serve as alternatives to the full formulation 
when faster time-integration is demanded. 


Remark. Though the reduced formulation with kf < 4 may also seem to be unconditionally stable in the 
above problem, the authors believe otherwise; we indeed were able to observe apparent free energy oscillation of 
O(10- 13 ) after steady state is roughly achieved around t — 4. This might be due to the numerical truncation 
error as the absolute residual error tolerance was set as 10 _1 ° throughout the analysis, but we further note 
that the full formulation with kf < 8 always produced solutions of decreasing free energy up to the precision 
of O(10 -16 ) even with the same level of error tolerance. The other observation of possible interest is that, the 
time-step size tends to put severer constraints on the solvability of the nonlinear system for reduced formulations 
as more high-order terms are ignored; it was more apparent for the reduced formulation with kf < 2 than for the 
reduced formulation with kf < 4. In the authors’ experience, however, those reduced formulations with kf < 4 
and kf < 2 still show better performance in this regard than the conventional Backward Euler method. 
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formulation 

t = 1 

t = 4 

ftp < 2 

12 

- 

ftp < 4 

16 

64 

ftp < 8 

23 

93 


Table 1: Time, in hour, required to compute solutions up to t = 1 and t = 4 on 32 s mesh with At = 5 x 10 -4 
using 2.60GHz Intel Xeon E5-2670 processors on 8 x 8 x 8 partitions for two reduced formulations, corresponding 
to ftp < 2 and ftp < 4, and the full formulation, corresponding to ftp < 8. 


ftp < 2 


ftp < 4 


ftp < 8 


* = 0.08 


* = 0.10 


* = 1.00 


* = 2.00 


* = 4.00 



Figure 9: Plots of e 2 and contour curves of c on deformed configurations at select times for two reduced formulations, 
one with ftp < 2 and the other one with ftp < 4, and the full formulation with ftp < 8. Solutions were computed 
on a 32 s mesh with At = 5 x 10 -4 . 




(a) (b) 

Figure 10: Plots of discrete total free energy IP” against time t for Q t £ [0,4] and |bj) f e [0.9,1.1] for two 
reduced formulations, one with ftp < 2 and the other one with ftp < 4, and the full formulation with ftp < 8. 
Solutions were computed on a 32 3 mesh with A* = 5x 10 -4 . 
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6 Conclusion 


We have developed a class of unconditionally stable, second-order accurate time-integration algorithms for 
nonlinear, mechano-chemical problems characterized by free energy functions that are non-convex in strain- 
composition space, and that must be stabilized by introducing spatial gradients of these fields. The associated 
phenomenology, which we term as mechano-chemical spinodal decomposition, includes the formation of mi- 
crostructural features and transient phenomena. Their resolution by first-order schemes such as the Backward 
Euler algorithm does not necessarily preserve the free energy decay that is a consequence of the second law of 
thermodynamics. 

The approach presented here has wide applicability in the design of stable, second-order schemes for coupled 
problems of mechanics and transport. Its use hinging on the Taylor expansion, can, in principle, be extended 
to any free energy function that is a multivariate polynomial of direct components, component gradients, and 
higher-order component gradients. Such functions are guaranteed to have finite Taylor-series, and therefore, 
analytic forms, which are required of the constitutive relations to preserve unconditional stability. While evalu¬ 
ation of the Taylor series expansions in the code comes at a cost, we have found that by exploiting symmetries 
inherent in the higher-order derivatives this cost can be substantially reduced. Furthermore, we have demon¬ 
strated that reduced formulations that truncate the higher-order terms in the Taylor series also can perform 
well for the initial and boundary value problems we have considered, even though unconditional stability can 
no longer be proven. Our analysis demonstrates, however, that second-order accuracy is maintained, regardless. 

To the best of our knowledge, this is the first treatment presenting stable and second-order schemes for 
systems that incorporate Toupin’s theory of gradient elasticity at finite strains. It has potential for extension 
to problems incorporating advection and reaction terms in the transport equation, as well as to systems that 
couple with the Allen-Cahn [T] treatment for evolution of non-conserved order parameters. It could thus cover 
a wide range of phase transformation phenomena involving solids as well as fluid phases in materials systems 
arising in battery, semi-conductor, polymer and structural applications. 
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